/*!
 ***************************************************************************
 *  \file nonuniform_meshes.cpp
 *  Using the StructuredMesh<> class with the NonUniform<> adaptor.
 ***************************************************************************
 *  \author Luis M. de la Cruz Salas [ Tue Dec 25 16:23:06 GMT 2010 ]
 ***************************************************************************
 *  This program presents examples for using the StructuredMesh<> class,
 *  which is intended to handle structured meshes in Cartesian coordinates.
 *  In these examples the class is specialized with the NonUniform<> adaptor 
 *  to create non-uniform meshes.
 *  This class have to be used in conjuntion with mesh type adaptors: 
 *  Uniform<> or NonUniform<> to create uniform or nonuniform meshes 
 *  respectively
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in
 *  according to your installation and then type the next commands:
 *  \verbatim 
    % make        
    % ./nonuniform_meshes \endverbatim
 *  The file \c input contains: 
 *  \verbatim 
    8  4 2  <--- Lenght in x, y and z axes respectively 
    11  9 6 <--- Number of nodes
    2       <--- Lambda factor \endverbatim
 *  The results of this program are written in the \c Data directory. 
 *
 *  In this example the function distribution is implemented by means of a
 *  functor, and is defined by the next function:
 *  \f[ 
 *  f(\xi) = { e^{\lambda\xi} - 1 \over e^{\lambda} - 1 } 
 *  \f]
 *  See the next figure:
 * \image html  f.png "Point distribution function for lambda = 2" width=10cm 
 * \image latex f.eps "Point distribution function for lambda = 2" width=10cm
 *
 *  The results of this example are shown in the next images. This figures 
 *  were generated with OpenDX, using the \c mesh.net script
 *  located in the \c Data directory.
 * \image html  mesh1D.png "1D Non-uniform mesh" width=10cm 
 * \image latex  mesh1D.eps "1D Non-uniform mesh" width=10cm 
 * \image html  mesh2D.png "2D Non-uniform mesh" width=10cm 
 * \image latex  mesh2D.eps "2D Non-uniform mesh" width=10cm 
 * \image html  mesh3D.png "3D Non-uniform mesh" width=10cm 
 * \image latex  mesh3D.eps "3D Non-uniform mesh" width=10cm 
 * \image html  mesh.png "2D Non-uniform mesh variable dist." width=10cm 
 * \image latex  mesh.eps "2D Non-uniform mesh variable dist." width=10cm .
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/StructuredMesh.hpp"
#include "Meshes/NonUniform.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Functor to define the point distribution.
class MyFunction {      
  double lambda; // distribution factor 
public:                                    
  MyFunction(double l) : lambda(l) { }   
  
  double operator() (double epsi) {      
    return ( (exp(lambda * epsi) - 1.0) / 
	     (exp(lambda       ) - 1.0) ); 
  }
};

// Scalar field in one dimension
typedef TunaArray<double,1>::huge ScalarField1D;
 
int main() 
{
    double length_x, length_y, length_z;
    int num_nodes_x, num_nodes_y, num_nodes_z;
    double lambda;

// Reading data from the input file
    std::ifstream input_cin ("input"); 
    input_cin >> length_x;
    input_cin >> length_y;
    input_cin >> length_z;
    input_cin >> num_nodes_x; 
    input_cin >> num_nodes_y; 
    input_cin >> num_nodes_z; 
    input_cin >> lambda;
    input_cin.close();

    StructuredMesh<NonUniform<double, 1> > mesh_1D(length_x, num_nodes_x);
    mesh_1D.calcPoints(X, MyFunction(lambda));
    mesh_1D.calcDeltas(X);
    mesh_1D.print();
    mesh_1D.writeToFile("./Data/mesh1D");

#ifdef WITH_GNUPLOT
    GNUplot plotter;
    plotter("p \"./Data/mesh1D\" w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
#endif    

    StructuredMesh<NonUniform<double, 2> > mesh_2D(length_x, num_nodes_x,
						   length_y, num_nodes_y);
    mesh_2D.calcPoints(X, MyFunction(lambda));
    mesh_2D.calcPoints(Y, MyFunction(lambda));
    mesh_2D.print();
    mesh_2D.writeToFile("./Data/mesh2D");

#ifdef WITH_GNUPLOT
    plotter("set view 0,0");
    plotter("splot \"./Data/mesh2D\" w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif

    StructuredMesh<NonUniform<double, 3> > mesh_3D(length_x, num_nodes_x,
						   length_y, num_nodes_y,
						   length_z, num_nodes_z);
    mesh_3D.calcPoints(X, MyFunction(lambda));
    mesh_3D.calcPoints(Y, MyFunction(lambda));
    mesh_3D.calcPoints(Z, MyFunction(lambda));
    mesh_3D.print();    
    mesh_3D.writeToFile("./Data/mesh3D");

#ifdef WITH_GNUPLOT
    plotter("set view 30,45");
    plotter("splot \"./Data/mesh3D\" using 1:2:3 w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif

//
//  2D Non-uniform mesh: It is possible make a different distribution on
//                       different parts of the mesh.
//
    StructuredMesh<NonUniform<double, 2> > mesh(length_x, num_nodes_x,
						length_y, num_nodes_y);
    int mid_nod_x = num_nodes_x / 2 + 1;
    int mid_nod_y = num_nodes_y / 2 + 1;
//
// x-axis distribution
//
    mesh.calcPoints(X, MyFunction(lambda),
		    0, length_x * 0.5,        // First half
		    mid_nod_x);            

    mesh.calcPoints(X, MyFunction(-lambda), 
		    length_x * 0.5, length_x, // Second half
		    mid_nod_x);            
					  
    mesh.calcDeltas(X);
//
// y-axis distribution
//
    mesh.calcPoints(Y, MyFunction(lambda),
		    0, length_y * 0.5,        // First half 
		    mid_nod_y);   
    
    mesh.calcPoints(Y, MyFunction(-lambda),
		    length_y * 0.5, length_y, // Second half
		    mid_nod_y);        
					  
    mesh.calcDeltas(Y);

    mesh.print();
    mesh.writeToFile("./Data/mesh");

#ifdef WITH_GNUPLOT
    plotter("set view 0,0");
    plotter("splot \"./Data/mesh\" w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif

    ScalarField1D x = mesh.getCoord(X);
    ScalarField1D y = mesh.getCoord(Y);
    ScalarField1D dxv = mesh.getDeltasVols(X);
    ScalarField1D dyv = mesh.getDeltasVols(Y);
    ScalarField1D dxf = mesh.getDeltasFace(X);    
    ScalarField1D dyf = mesh.getDeltasFace(Y);

    std::cout << "\n >---> Some info: \n";
    std::cout << "\n x = " << x << "\n";
    std::cout << "\n y = " << y << "\n";
    std::cout << "\n dvols x = " << dxv << "\n";
    std::cout << "\n dvols y = " << dyv << "\n";
    std::cout << "\n dface x = " << dxf << "\n";
    std::cout << "\n dface y = " << dyf << "\n";
    
    return 0;
}
